/*
 * Edge Case Handling for GPU Correctness Verification
 * Handles zero arrays, NaN, Infinity, and denormal float values
 */

#include <cuda_runtime.h>
#include <math.h>
#include <cstdio.h>
#include <cstdlib.h>
#include <string>
#include <vector>

// Edge case types
enum EdgeCaseType {
    EDGE_ZERO = 0,
    EDGE_NAN = 1,
    EDGE_INFINITY = 2,
    EDGE_DENORMAL = 3,
    EDGE_LARGE = 4,
    EDGE_SMALL = 5,
    EDGE_MIXED = 6
};

// Edge case handling configuration
struct EdgeCaseConfig {
    EdgeCaseType type = EDGE_ZERO;
    float large_value = 1e6f;
    float small_value = 1e-6f;
    float denormal_threshold = 1e-38f;
    bool check_special_values = true;
    bool generate_report = true;
    int nan_frequency = 100;  // Every Nth element
    int inf_frequency = 50;   // Every Nth element
};

// Edge case validation result
struct EdgeCaseResult {
    bool passed = false;
    int total_elements = 0;
    int zero_count = 0;
    int nan_count = 0;
    int inf_count = 0;
    int denormal_count = 0;
    int large_count = 0;
    int small_count = 0;
    std::vector<int> error_indices;
    std::vector<float> error_values;
    
    void print_summary() const {
        printf("Edge Case Validation Summary:\n");
        printf("  Total elements: %d\n", total_elements);
        printf("  Zero values: %d\n", zero_count);
        printf("  NaN values: %d\n", nan_count);
        printf("  Infinity values: %d\n", inf_count);
        printf("  Denormal values: %d\n", denormal_count);
        printf("  Large values: %d\n", large_count);
        printf("  Small values: %d\n", small_count);
        printf("  Errors: %d\n", (int)error_indices.size());
    }
};

// GPU kernel for edge case handling
__global__ void edge_case_kernel(
    const float* __restrict__ a,
    const float* __restrict__ b,
    float* __restrict__ c,
    int n,
    int edge_case_type,
    int* __restrict__ status_flags,
    float* __restrict__ debug_info,
    int debug_offset
) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    
    if (idx < n) {
        float a_val = a[idx];
        float b_val = b[idx];
        float result;
        
        switch (edge_case_type) {
            case EDGE_ZERO:
                result = 0.0f + 0.0f;
                break;
                
            case EDGE_NAN:
                // Handle NaN propagation
                if (isnan(a_val) || isnan(b_val)) {
                    result = NAN;
                } else {
                    result = a_val + b_val;
                }
                break;
                
            case EDGE_INFINITY:
                // Handle infinity propagation
                if (isinf(a_val) || isinf(b_val)) {
                    if (isinf(a_val) && isinf(b_val) && signbit(a_val) != signbit(b_val)) {
                        result = NAN;  // inf - inf = NaN
                    } else {
                        result = a_val + b_val;
                    }
                } else {
                    result = a_val + b_val;
                }
                break;
                
            case EDGE_DENORMAL:
                // Handle denormals by flushing to zero
                a_val = (fabsf(a_val) < 1e-38f) ? 0.0f : a_val;
                b_val = (fabsf(b_val) < 1e-38f) ? 0.0f : b_val;
                result = a_val + b_val;
                break;
                
            case EDGE_LARGE:
                // Handle large values with overflow checking
                result = a_val + b_val;
                if (isinf(result)) {
                    result = INFINITY;
                }
                break;
                
            case EDGE_SMALL:
                // Handle small values with underflow checking
                result = a_val + b_val;
                if (fabsf(result) < 1e-38f && result != 0.0f) {
                    result = 0.0f;  // Flush to zero
                }
                break;
                
            case EDGE_MIXED:
                // Mixed edge cases based on index
                if (idx % 100 == 0) {
                    result = NAN;
                } else if (idx % 50 == 0) {
                    result = INFINITY;
                } else if (idx % 25 == 0) {
                    result = 0.0f;
                } else {
                    result = a_val + b_val;
                }
                break;
                
            default:
                result = a_val + b_val;
        }
        
        c[idx] = result;
        
        // Store status information
        if (status_flags) {
            int flag = 0;
            if (isnan(result)) flag |= 1;
            if (isinf(result)) flag |= 2;
            if (fabsf(result) < 1e-38f && result != 0.0f) flag |= 4;
            if (result == 0.0f) flag |= 8;
            status_flags[idx] = flag;
        }
        
        // Store debug information
        if (debug_info && idx < debug_offset) {
            debug_info[idx * 8] = a_val;
            debug_info[idx * 8 + 1] = b_val;
            debug_info[idx * 8 + 2] = result;
            debug_info[idx * 8 + 3] = (float)idx;
            debug_info[idx * 8 + 4] = isnan(result) ? 1.0f : 0.0f;
            debug_info[idx * 8 + 5] = isinf(result) ? 1.0f : 0.0f;
            debug_info[idx * 8 + 6] = (fabsf(result) < 1e-38f && result != 0.0f) ? 1.0f : 0.0f;
            debug_info[idx * 8 + 7] = (result == 0.0f) ? 1.0f : 0.0f;
        }
    }
}

// Special value validation kernel
__global__ void special_value_validation_kernel(
    const float* __restrict__ computed,
    const float* __restrict__ expected,
    int* __restrict__ error_flags,
    float* __restrict__ error_values,
    int n,
    float tolerance
) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    
    if (idx < n) {
        float comp = computed[idx];
        float exp = expected[idx];
        
        // Special value handling
        bool is_comp_nan = isnan(comp);
        bool is_exp_nan = isnan(exp);
        bool is_comp_inf = isinf(comp);
        bool is_exp_inf = isinf(exp);
        
        int error = 0;
        
        if (is_comp_nan != is_exp_nan) {
            error = 1;  // NaN mismatch
        } else if (is_comp_inf != is_exp_inf) {
            error = 2;  // Infinity mismatch
        } else if (!is_comp_nan && !is_comp_inf && !is_exp_nan && !is_exp_inf) {
            float diff = fabsf(comp - exp);
            if (diff > tolerance) {
                error = 3;  // Value mismatch
            }
        }
        
        if (error_flags) {
            error_flags[idx] = error;
        }
        
        if (error_values && error != 0) {
            error_values[idx * 2] = comp;
            error_values[idx * 2 + 1] = exp;
        }
    }
}

// Edge case generator class
class EdgeCaseGenerator {
private:
    EdgeCaseConfig config;
    
public:
    EdgeCaseGenerator(const EdgeCaseConfig& cfg = EdgeCaseConfig()) : config(cfg) {}
    
    // Generate edge case data
    void generate_edge_case_data(float* a, float* b, int size, 
                               EdgeCaseType type = EDGE_MIXED) {
        switch (type) {
            case EDGE_ZERO:
                generate_zero_arrays(a, b, size);
                break;
            case EDGE_NAN:
                generate_nan_arrays(a, b, size);
                break;
            case EDGE_INFINITY:
                generate_infinity_arrays(a, b, size);
                break;
            case EDGE_DENORMAL:
                generate_denormal_arrays(a, b, size);
                break;
            case EDGE_LARGE:
                generate_large_arrays(a, b, size);
                break;
            case EDGE_SMALL:
                generate_small_arrays(a, b, size);
                break;
            case EDGE_MIXED:
                generate_mixed_arrays(a, b, size);
                break;
        }
    }
    
    void generate_zero_arrays(float* a, float* b, int size) {
        for (int i = 0; i < size; i++) {
            a[i] = 0.0f;
            b[i] = 0.0f;
        }
    }
    
    void generate_nan_arrays(float* a, float* b, int size) {
        for (int i = 0; i < size; i++) {
            if (i % config.nan_frequency == 0) {
                a[i] = NAN;
                b[i] = 1.0f;
            } else {
                a[i] = (float)(i % 100);
                b[i] = (float)(i % 100);
            }
        }
    }
    
    void generate_infinity_arrays(float* a, float* b, int size) {
        for (int i = 0; i < size; i++) {
            if (i % config.inf_frequency == 0) {
                a[i] = (i % 2 == 0) ? INFINITY : -INFINITY;
                b[i] = 1.0f;
            } else {
                a[i] = (float)(i % 1000);
                b[i] = (float)(i % 1000);
            }
        }
    }
    
    void generate_denormal_arrays(float* a, float* b, int size) {
        for (int i = 0; i < size; i++) {
            if (i % 200 == 0) {
                a[i] = 1e-39f;  // Denormal
                b[i] = 1e-39f;
            } else {
                a[i] = (float)(i % 1000);
                b[i] = (float)(i % 1000);
            }
        }
    }
    
    void generate_large_arrays(float* a, float* b, int size) {
        for (int i = 0; i < size; i++) {
            a[i] = (i % 2 == 0) ? config.large_value : -config.large_value;
            b[i] = (i % 2 == 0) ? config.large_value : -config.large_value;
        }
    }
    
    void generate_small_arrays(float* a, float* b, int size) {
        for (int i = 0; i < size; i++) {
            a[i] = (i % 2 == 0) ? config.small_value : -config.small_value;
            b[i] = (i % 2 == 0) ? config.small_value : -config.small_value;
        }
    }
    
    void generate_mixed_arrays(float* a, float* b, int size) {
        for (int i = 0; i < size; i++) {
            int pattern = i % 1000;
            
            if (pattern < 100) {
                a[i] = 0.0f;
                b[i] = 0.0f;
            } else if (pattern < 200) {
                a[i] = NAN;
                b[i] = 1.0f;
            } else if (pattern < 300) {
                a[i] = (pattern < 250) ? INFINITY : -INFINITY;
                b[i] = 1.0f;
            } else if (pattern < 400) {
                a[i] = 1e-39f;
                b[i] = 1e-39f;
            } else if (pattern < 500) {
                a[i] = config.large_value;
                b[i] = config.large_value;
            } else if (pattern < 600) {
                a[i] = config.small_value;
                b[i] = config.small_value;
            } else {
                a[i] = (float)(pattern - 600);
                b[i] = (float)(pattern - 600);
            }
        }
    }
    
    // Validate edge case results
    EdgeCaseResult validate_edge_cases(const float* computed, const float* expected, int size,
                                     float tolerance = 1e-6f) {
        EdgeCaseResult result;
        result.total_elements = size;
        
        for (int i = 0; i < size; i++) {
            float comp = computed[i];
            float exp = expected[i];
            
            // Count special values
            if (comp == 0.0f) result.zero_count++;
            if (isnan(comp)) result.nan_count++;
            if (isinf(comp)) result.inf_count++;
            if (fabsf(comp) < config.denormal_threshold && comp != 0.0f) 
                result.denormal_count++;
            if (fabsf(comp) > 1e6f) result.large_count++;
            if (fabsf(comp) < 1e-6f && comp != 0.0f) result.small_count++;
            
            // Validate against expected
            bool is_nan_comp = isnan(comp);
            bool is_nan_exp = isnan(exp);
            bool is_inf_comp = isinf(comp);
            bool is_inf_exp = isinf(exp);
            
            bool valid = true;
            
            if (is_nan_comp != is_nan_exp) {
                valid = false;
            } else if (is_inf_comp != is_inf_exp) {
                valid = false;
            } else if (!is_nan_comp && !is_inf_comp) {
                float diff = fabsf(comp - exp);
                if (diff > tolerance) valid = false;
            }
            
            if (!valid) {
                result.error_indices.push_back(i);
                result.error_values.push_back(comp);
            }
        }
        
        result.error_count = result.error_indices.size();
        result.passed = (result.error_count == 0);
        
        if (config.generate_report) {
            result.print_summary();
        }
        
        return result;
    }
    
    // Launch edge case kernel
    void launch_edge_case_kernel(
        const float* d_a, const float* d_b, float* d_c, int n,
        EdgeCaseType type, int* d_status_flags, float* d_debug_info,
        int threads_per_block = 256
    ) {
        int blocks = (n + threads_per_block - 1) / threads_per_block;
        edge_case_kernel<<<blocks, threads_per_block>>>(
            d_a, d_b, d_c, n, type, d_status_flags, d_debug_info, 1000
        );
    }
};

// Edge case test suite
class EdgeCaseTestSuite {
public:
    static bool run_all_tests(int array_size, float tolerance = 1e-6f) {
        bool all_passed = true;
        
        EdgeCaseConfig config;
        EdgeCaseGenerator generator(config);
        
        std::vector<float> a(array_size);
        std::vector<float> b(array_size);
        std::vector<float> expected(array_size);
        std::vector<float> computed(array_size);
        
        // Test zero arrays
        generator.generate_edge_case_data(a.data(), b.data(), array_size, EDGE_ZERO);
        for (int i = 0; i < array_size; i++) {
            expected[i] = a[i] + b[i];
        }
        
        if (!test_case("Zero Arrays", a.data(), b.data(), expected.data(), 
                      computed.data(), array_size, tolerance)) {
            all_passed = false;
        }
        
        // Test NaN arrays
        generator.generate_edge_case_data(a.data(), b.data(), array_size, EDGE_NAN);
        for (int i = 0; i < array_size; i++) {
            if (isnan(a[i]) || isnan(b[i])) {
                expected[i] = NAN;
            } else {
                expected[i] = a[i] + b[i];
            }
        }
        
        if (!test_case("NaN Arrays", a.data(), b.data(), expected.data(),
                      computed.data(), array_size, tolerance)) {
            all_passed = false;
        }
        
        // Test infinity arrays
        generator.generate_edge_case_data(a.data(), b.data(), array_size, EDGE_INFINITY);
        for (int i = 0; i < array_size; i++) {
            if (isinf(a[i]) || isinf(b[i])) {
                if (isinf(a[i]) && isinf(b[i]) && signbit(a[i]) != signbit(b[i])) {
                    expected[i] = NAN;
                } else {
                    expected[i] = a[i] + b[i];
                }
            } else {
                expected[i] = a[i] + b[i];
            }
        }
        
        if (!test_case("Infinity Arrays", a.data(), b.data(), expected.data(),
                      computed.data(), array_size, tolerance)) {
            all_passed = false;
        }
        
        return all_passed;
    }
    
private:
    static bool test_case(const char* name, const float* a, const float* b,
                         const float* expected, float* computed, int size, float tolerance) {
        printf("Testing %s...\n", name);
        
        // Simple CPU computation for testing
        for (int i = 0; i < size; i++) {
            computed[i] = a[i] + b[i];
        }
        
        EdgeCaseConfig config;
        EdgeCaseGenerator generator(config);
        
        EdgeCaseResult result = generator.validate_edge_cases(computed, expected, size, tolerance);
        
        printf("  %s: %s (%d errors)\n", name, 
               result.passed ? "PASSED" : "FAILED", result.error_count);
        
        return result.passed;
    }
};

// C-style interface functions
extern "C" {
    EdgeCaseGenerator* create_edge_case_generator(int case_type, float large_value, 
                                                  float small_value, float denormal_threshold) {
        EdgeCaseConfig config;
        config.type = static_cast<EdgeCaseType>(case_type);
        config.large_value = large_value;
        config.small_value = small_value;
        config.denormal_threshold = denormal_threshold;
        return new EdgeCaseGenerator(config);
    }
    
    void destroy_edge_case_generator(EdgeCaseGenerator* generator) {
        delete generator;
    }
    
    void generate_edge_case_data(float* a, float* b, int size, int case_type) {
        EdgeCaseConfig config;
        EdgeCaseGenerator generator(config);
        generator.generate_edge_case_data(a, b, size, static_cast<EdgeCaseType>(case_type));
    }
    
    int validate_edge_cases(const float* computed, const float* expected, int size,
                           float tolerance, int* error_count) {
        EdgeCaseConfig config;
        EdgeCaseGenerator generator(config);
        EdgeCaseResult result = generator.validate_edge_cases(computed, expected, size, tolerance);
        
        if (error_count) *error_count = result.error_count;
        return result.passed ? 1 : 0;
    }
    
    int run_edge_case_tests(int array_size, float tolerance) {
        return EdgeCaseTestSuite::run_all_tests(array_size, tolerance) ? 1 : 0;
    }
}